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Abstract 

Starting from a travelling wave ansatz we show successively that the shape 
of a nonlinear excitation generally depends also on the 1 st , 2 nd , . . . time deriva- 
tive of the position X of the excitation. From the Hamilton equations we derive 
a hierarchy of equations of motion for X. The type of the excitation determines 
on which levels the hierarchy can be truncated consistently: "Gyrotropic" exci- 
tations are governed by odd-order equations, non-gyrotropic ones by even-order 
equations. Examples for the latter case are kinks in 1-dimensional models and 
planar vortices of the 2-d anisotropic (easy-plane) Heisenberg model. The non- 
planar vortices of this model are the simplest gyrotropic example. For this case 
we solve the Hamilton equations for a finite system with one vortex and free 
boundary conditions and calculate the parameters of the 3 rd -order equation 
of motion. This equation yields trajectories which are a superposition of two 
cycloids with different frequencies, which is in full agreement with computer 
simulations of the full many-spin model. Finally we demonstrate that the 
additional effects from the 5 th -order equation are negligible. 

1 Introduction 

Nonlinear coherent excitations, such as solitons or solitary waves, usually have some 
particle-like properties. E.g., the equation of motion of their 'center-of-mass' is New- 
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tonian, at least in a first approximation. For the case of zero force the excitations 
can move at constant velocity due to their inertial mass. However, there are also 
other excitations, which do not behave like classical particles, i.e., Galilei's law is not 
valid. For spin systems, good examples are certain magnetic domains and non-planar 
vortices in two- or three-dimensional models. For this type of nonlinear collective ex- 
citations Thiele |], derived a l st -order equation of motion for the collective variable 
X(t) 

XxG = F, (1) 

which we will refer to as the Thiele Eq. in the following. It is valid only for steady- 
state motion because it was derived from the Landau-Lifshitz equation for spins 
assuming a rigid shape of the excitations. 

F is either an external force or the force due to interactions with other excitations. 
V x G is called a gyrocoupling force. It is formally equivalent to the Lorentz force. 
However G is not an external magnetic field but an intrinsic quantity, produced by 
the excitation itself and carried along with it. G is called the gyrocoupling vector, or 
for short gyrovector. We call excitations with \G\ = G ^ gyrotropic. For 1-d spin 
models G is always zero. For the vortices of the 2-d anisotropic Heisenberg model 

G = 27cqpe z (2) 

as was first calculated by Huber |§. Here e z is the unit vector perpendicular to the 
xy-plane in which the magnetic ions are situated, q = ±1, ±2, ... is the vorticity, p 
is a second topological charge which is defined as the value of the spin component 
S z at the vortex center in the continuum limit. 

The 2-d anisotropic Heisenberg model with XF-symmetry is defined by the spin 
Hamiltonian 

H = —J 2^ Pm'-'n + SmS n + (1 — "J^m^nl (^) 

<m,n> 

with < S < 1. Here <m,n> labels the nearest-neighbor sites of a square lattice. 
We treat the spin S as a classical vector and set S = J = 1 . 

Due to the anisotropy S the spins prefer to be oriented in the xy-p\a.ne which is 
therefore often called easy plane. For strong anisotropy (6 C < 5 < 1, with S c ~ 0.28) 
only planar vortices are stable @, ||. In the static case their spins are lying completely 
in the plane, while dynamically small S^-components develop. But at the vortex 
center S z is always zero, thus p = and G = due to (|2|). Therefore the Thiele Eq. 
is not valid here. 

In contrast to this case, for weak anisotropy (0 < 5 < 5 C ) only non-planar vortices 
are stable 0. They exhibit a localized structure of the S z - components around the 
vortex center, at which S z = ±1. Although this structure changes due to motion, the 
value at the center remains the same (section [j|). Thus p = ±1 determines to which 
side the out-of-plane structure of the vortex is oriented and is therefore be termed 
polarization. The Thiele Eq. is valid for steady-state motion when the vortex shape 



is rigid (in the moving frame). This includes the case of a constant rotation on a 
circle. 

In 1994 Wysin et al. dropped the rigid-shape assumption by allowing the 

vortex shape to depend on the velocity X(t) and derived a generalized Thiele equation 

MX + X xG = F (4) 

As the velocity dependent parts of the vortex structure decay like 1/r with the 
distance r from the vortex center J3J, the mass was predicted to be 

M~^ln(L/a ) (5) 

where L and a are upper and lower cut-offs in the order of the system size and lattice 
constant, respectively. The same mass also appears in the kinetic energy MX 2 /2 of 
a vortex j4|. 

According to Eq. Q) the trajectory X(t) of the vortex center is formally the same 
as that of an electric charge e in a plane with a perpendicular magnetic field B and 
an in-plane electric force F, namely a cycloid with frequency uj = G/M, cf. the 
cyclotron frequency eB/(Mc), where c is the speed of light. However, a test of this 
prediction by computer simulations turned out to be rather difficult: For the model 
(0) the use of an external force creates additional nonlinear excitations; e.g., an in- 
plane magnetic field creates a double domain wall which connects the vortex with a 
boundary [0. Therefore 2- vortex simulations were performed where each vortex is 
driven by the force between them || [|. Large square systems with free boundaries 
were used and the trajectories were chosen such that the vortices stayed far away 
from each other and from the boundaries. 

A 2- vortex theory was developed §J which explains one very important qualitative 
feature of the simulations: There are four main scenarios in vortex dynamics, however 
cycloidal oscillations around a mean trajectory X°(t) can be observed only for two of 
them, namely where either a vortex and an antivortex (i.e. q\ = — q%) rotate around 
each other, or where two vortices with equal vorticity perform a parallel translation 
with constant average speed. There are no oscillations for the other 2 scenarios of 
'vortex- vortex rotation' and 'vortex-antivortex translation'. 

However, the frequency u of the observed cycloidal oscillations yielded a massQ 
M = G/uj = 2-k/uj which was much larger than predicted by Eq. (|5J), see Table lb 
of Ref . [f| . Later a second severe discrepancy was discovered after having improved 
the simulations for the 2- vortex rotation |TD] : Better initial conditions were used and 



instead of a square a circular system was used, with symmetric initial positions. Here 
for each vortex the mean trajectory X°(t) is a circle; the spectrum of the oscillations 
can be measured much more accurately and shows very clearly two frequencies cj 1>2 

^^More precisely, for a 2-vortex system an effective mass tensor was defined, its eigenvalues were 
both much larger than predicted. 



with about the same strength, instead of the one frequency uj = G/M predicted by 
Eq. (§. 

Because of the above two discrepancies between theory and simulations we now 
derive a new theory which amounts to a hierarchy of equations of motion for a 
nonlinear coherent excitation in a system with an arbitrary Hamiltonian (section 0). 
The derivation is completely general, but for simplicity we consider a Hamiltonian 
which is a functional of only one field and its canonical momentum; this is the case 
of our model (|^) in the continuum limit. 

The level n in this hierarchy is defined as the order n of the highest time derivative 
which appears in the equation of motion. A classification of the excitations deter- 
mines on which levels the hierarchy can be truncated consistently: The dynamics of 
gyrotropic excitations is described by odd-order equations, the simplest example is 
the non-planar vortex of model (||). Non-gyrotropic excitations (e.g. kinks in 1-d 
models, planar vortices in model (|^)) are governed by even-order equations of motion. 

In order to obtain the above classification one must know at least the order of 
magnitude of the parameters in the equations of motion. For the calculation of the 
parameters it is necessary to solve the Hamilton equations for a system with one 
excitation, choosing appropriate boundary conditions. We do this for our model ([|) 
in section |] using a finite system with free boundaries. 

After the calculation of the parameters the equation of motion can be solved. For 
the non-planar vortices we solve a 3 rd -order equation, which in fact yields the observed 
frequencies Ui t 2 (the trivial solution u = yields the mean trajectory X°(t)), section 
|j. We compare with new simulations where a single vortex on a circle is driven by its 
image vortex (section |5|). We only note that the 2- vortex simulations yield a similar 
spectrum but cannot be compared quantitatively with this 1-vortex theory. We also 
show that the 5 th -order equation, which appears on the next consistent level of the 
hierarchy, predicts two additional frequencies which can also be observed, but they 
are very weak. Therefore all levels higher than three can safely be neglected in this 
case. 

In section |5| we discuss our results and in the appendix we describe the numerical 
part of this work in more detail. 

2 Hierarchy of equations of motion 

We consider an arbitrary classical Hamiltonian which is a functional of a field </>(f, t) 
and its canonical momentum m(r,t). We introduce a collective variable X(t) for 
the position of a nonlinear coherent excitation. We make a travelling wave ansatz 
4>(f,t) = (j)(f — X(t)), m(f,t) = m(r — X(t)) where, strictly speaking, the functions 
on the r.h.s. should bear an index, which is omitted here for simplicity. We insert 

X i ' rh = l^ X i ' ( 6 ) 



dXj J ' dX d 



where a summation over j — 1, 2, 3 is implied, into the Hamilton equations 

= ^ ' m = ~w (7) 

Since the r.h.s. of these equations generally contain m and <fi it is clear that m and 

4> depend not only on X but also on X. Thus the shape of the collective excitations 
generally depends on the velocity and we take this into account by the improved 

ansatz <j)(f,t) = <fi(r — X,X), m(f,t) = m[f — X,X). However, if we now insert 

and m into the eqs. (|7|) we see an additional dependence on X, and so on. Obviously 
we must truncate somewhere. We will see below that a truncation at the derivative 
X 1 -™™ 1 ) yields an n th -order equation of motion. For gyrotropic excitations only odd- 
order equations will turn out to be consistent (section £|). As an example we derive 
here the 3 rd -order equation by the ansatz 

<f>(r,t) = <f>(r-X(t)j(t),X(t)) (8a) 

m(f,t) = m(r-X(t),X(t),X(t)). (8b) 

It is unusual that a collective excitation should depend on the acceleration, but 
it will be seen later that this dependence yields one of the dominant terms in the 
equation of motion. This equation is now easily derived by using the technique of 
Wysin et al. ||. Our ansatz is inserted into the Hamilton equations (0) and yields 

d(j) ■ d<f) ~ d<f) ... 5H 

dXj J dXj 3 dXj 3 5m 

dm ■ dm ■■ dm ■■■ 5H , , , 

a^ x ' + ax- x ' + W, x ' = ~W (9b) 

where H depends on X, X and X via (j) and m. Multiplying these equations by ^^ 
and -g^r, resp., subtracting the equations from each other, and integrating over f we 
obtain 

A X +MX + GX = F (10) 

Here 

with the Hamiltonian density 7i is the static force, it is either an external force or 
the force due to the interactions with other nonlinear excitations. 

3 J [dXidXj dXjdXij K ' 

is equivalent to the gyrocoupling tensor of Thiele Jl|, §]. Because of the antisymmetry 
of G the term GX can be written as —G x X. The gyrovector G is orthogonal to 



the plane denned by X and X and has already been discussed in the introduction. 

The terms GX and F constitute the Thiele Eq. valid only for steady-state motion. 
We note that the above method is much shorter than the original one by Thiele, who 
started directly from the Landau-Lifshitz Eq. (|5^|). 

f ,s (96 dm 86 dm } 
13 J [dXidXj dXjdXij v ; 

is referred to as the mass tensor, and 

f * (96 dm 86 dm 1 

An = I drr { — - — r. i-- — } (14) 

as the 3 rd -order gyrotensor. In order to see which properties they have it is necessary 
to calculate the functions 6 and m in the ansatz (|j) which describe the shape of the 
excitations. We do this for our model (Bf). 



3 Hamilton equations for a finite system with one 
vortex and free boundaries 

We consider the 2-d anisotropic Heisenberg model (H) and introduce 

<K?,t) = tan_1 ||r§ ( 15a ) 

m(r,t) = S z (f,t), (15b) 

where the site index at the spin components has been replaced by the variable r = 
{xxiX%\. The Hamiltonian is 



(Vm) 



H=- fd 2 rl(l- m 2 ) ( V0) 2 + 5[Am 2 - ( Vm) 2 ] + -^ 
and the Hamilton equations read 



m? 



(16) 



^ m[4Hwf ,(J_. j)im .^ (17a) 

— = (l-m 2 )A6-2mVmV6. (17b) 

For < 5 < 5 C (with S c ~ 0.28 for a square lattice) only non-planar vortices are 
stable and they have the following static structure Q 

6 = gtan" 1 - 2 --^^ (18) 
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for r' <C ry 



with 



r-X. 



(19) 
(20) 

(21) 



The polarization isp = ±1 and only vorticities q — ±1 are considered. The constants 
ai and a 2 can be determined by matching the two solutions ( |I9"D and fl20| ) at some 
intermediate distance, e. g. at r « 7y. The characteristic length 



?V 




is interpreted as the radius of the vortex core. The constant (fo in 



From the static structure one obtains for (|T2"D 



G 



i.i 



Geij 



(22) 
stays arbitrary. 

(23) 



with G = 2irqp. eij is the 2-d completely antisymmetric tensor. 
For a slowly moving vortex we expand in a perturbation series 

4> = cf> (r') + 4>i(r'J) + Mr'J) 
m = mo(r')+mi(f",X)+m2(r',X) 



(24a) 
(24b) 



and assume that the first and second order terms depend linearly on Xj and Xj, 
resp. Slow motion means here that the velocity is much smaller than the spin-wave 
velocity which is 2y5 in our dimensionless units 0. 

In the next section it will be seen that the main contributions to the integrals for 
M and A stem from the region far away from the vortex center, i.e. from r' 3> ry. 
Particularly the dependence on the system size naturally comes entirely from this 
region. Therefore it is sufficient to solve the Hamilton equations only for this region. 



We assume that here only the terms 45m and A0 on the r.h.s. of (|17a|) and (|17b|) , 
resp., are important for the dynamic parts of and m. This assumption can be 
justified a posteriori. 

In O(Xj) we then obtain for r' ^> ry 



<90o 
dX j 



Xj = 4Smi 



dm 
dX~ 



The first equation gives 



mi 



48(r' 



— ( x 2 X\ — x x X 2 



Xj = A0i. 



A5r' ' '' 



(25) 



(26) 



where (p 1 is the angle between r ' and X. For the second equation in (|25| ) the inho- 
mogenity vanishes at large distances. We therefore have the solution 

0! = pcx'jXj (27) 

for large r' which increases linearly towards the free boundaries, in contrast to the 
asymptotic solution of ref. [|J where an infinite system with decaying boundary con- 
ditions was considered. The constant c can be determined by matching the solution 
([27|) to the small-r' solution of ref. |4j] . 

In 0(Xj) we obtain by the above methods 

m 2 = pbx'jXj (28) 

with b = c/(45), and 

fa = -jL]nr , (x , 2 X 1 -x' 1 X 2 ). (29) 

Here m 2 ~ cosx' and 2 ~ sinx', where x' is the angle between r' and X. 

The results of this section can be tested by looking at snapshots of the orientation 
of the spins in our computer simulations for the discrete system (the simulations are 
discussed in section |5|). The dynamic parts of are very difficult to observe because 
(f) must be substracted first. But this depends very sensitively on the vortex position, 
which is known only within a certain accuracy. However, mo vanishes exponentially 
for r' 3> ry. Here the dynamic parts of m can be observed directly. Though m 2 
increases linearly with r' while m x decays, m 2 generally cannot be distinguished 

clearly from mi because generally \X\ -C \X\. In order to observe m 2 nevertheless, 
we have selected specific points of the vortex trajectories in section ^[ At the turning 
points the acceleration has a maximum while the velocity is small. Here the predicted 
cos x'-dependence of m 2 and its linear increase towards the boundaries are seen very 
clearly (Figs. [I] and |3|). In the middle between two turning points the trajectory 
is nearly straight and the acceleration is small. Here m 2 is in fact barely visible 
(Fig. 0), instead several humps can be observed which are probably produced by mi. 
By the evaluation of contour plots of the above snapshots we estimate b ~ 2.5 for 
the parameter in (^). This parameter is the only one we need for the calculation of 
the integrals in the following section. 

4 Mass, third-order gyrotensor and solution of the 
equation of motion 

Since we consider very slow motion we need to calculate only the rest mass, i.e. in 
Eq. flTB| ) the derivations with respect to Xi are applied only to the static parts of <p 
and m. As </>i and mi depend linearly on the velocity, 

13 J \ dX % dXj dXj dXif K J 



is the rest mass. In the same way we obtain the constant part of (|TJ[) by 

a _ fj2 r \d<f> dm 2 <90 2 <9m o l 

An analytic calculation is possible if we choose a circular system (radius L) and 
consider a vortex with its center A at the circle center, which is chosen as origin. 
We divide the integration region into an inner part < r < a c and an outer part 
& c < r < L, where we choose a c 3> ry. The inner regions yield L-independent 
contributions, while the contributions from the outer region turn out to increase 



with L and thus dominate for large L. As mo decays exponentially according to (|19D , 
the second terms in (|30|) and fl3~l|) give no contribution for the outer region. Therefore 
we need to calculate only the first terms and obtain with (fL8|), (f26[) and ( ^8|) 

7T L 

Mu = MS- , M = — In — + const. (32a) 

Aij = Aeij , A = — {l 2 - a 2 ^j + const. (32b) 

where (^ is the 2-d unit matrix. The constants are the above L-independent con- 
tributions from the inner region. They depend on a c , but together with the other 
a c -dependent terms they must add up to a c -independent constants M and A be- 
cause the final result 

M = i-AnL + M , A = —L 2 + A (33) 

4o 4 

must not depend on how the integration region is divided. 

If the vortex center is not at the circle center the situation is less symmetric, but 
M is still diagonal (in a system with axes parallel and orthogonal to X) and the An 
are still zero. The non-vanishing matrix elements must be calculated by numerical 
integration. However, in our simulations the mean vortex trajectory is a circle with 
radius Rq ^ L around the circle center. Therefore the differences between M u and 
M 2 2 and those between A 12 and — A 2 \ are small and will be neglected in the following. 

For the solution of the equation of motion ([[(]) we consider a small displacement 
from a mean trajectory X°(t) 

X(t)=X°(t)+x(t). (34) 

We consider a situation where the force is always pointing in the Ax-direction and 
expand to 1 st order around A° = R 

F = F + F" oXl . (35) 

The mean trajectory is a parallel to the A 2 -axis, namely A° = R , X® = Vot, where 
Vq = Fq/G is positive for a vortex with qp = 1. 



For the displacement x we get two linear equations 



which are solved by 



A x 2 +Mxi + Gx 2 = FqX X (36a) 

- A x\ +Mx 2 -Gxi = (36b) 



x\ = ai cos u>it + a,2 cos u> 2 t (37a) 

x 2 = h\ sin uj\t + b 2 sin u 2 t (37b) 



with 



2 2GA + M 2 (AGA + M 2 )M 2 MF^ 
^ = 2A 2 T V AA* + ^4^- (38) 



The amplitude ratios K a = b a /a a are obtained from 

(G-Awl)w a K a = Mu 2 a + F^ (39a) 

G-Au 2 a = Mu a n a (39b) 

which yields 

n a = ±sjl + FM(Mul), a = 1,2. (40) 

In order to facilitate the discussion we set Fq = (our simulations are made anyway 
for the case of small Fq and Fq. Then we get 



G fMV M 
^1,2 = \ 



\ 2 M 
A ' \2A) T 2l (41) 

with Ki 5 2 = ±1. As A ~ L 2 while M ~ InL for a large system, the leading In- 
dependence of ui a is 1/L . 

The frequencies 0*1,2 form a narrow doublet. The mean frequency depends on A, 
but not on M 

u c = v^ZJs" = y/G/A ~ 1/L. (42) 

Contrary to this, the splitting of the doublet 

InL 

Au = uj 2 - u x = MIA — (43) 

is proportional to the mass. 

We now discuss the size dependence of the different terms in the equation of 
motion (|To|) . As every time derivative contributes a factor u a ~ 1/L, we see that 
A Xj~ 1/L and MXi ~ In L/L 2 for a large system. Therefore A Xi cannot be 
neglected when MXi is retained. This is the reason for the discrepancies resulting 
from the 2 nd -order equation (|j), see the introduction. However, the neglection of 
both AXi and MXi represents a consistent approximation, namely the Thiele Eq. 

(ED- 
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We have just seen that the levels n — 1 and n = 3 represent consistent approx- 
imations in a hierarchy of equations of motion. Let us now include two additional 
time derivatives in our ansatz (^) which yields two additional orders in the equation 
of motion. Using the methods of sections |3] and |] one can show that the parameters 
in the 5 th - and 4 th -order terms scale like L 4 and L 2 InL, resp.. The same arguments 
as above show that n = 5 is the next consistent level. Here we get 4 frequencies u a 
(besides the trivial solution u = 0). Neglecting for the moment the small 2 nd - and 
4 th -order terms, we obtain two two-fold degenerate frequencies. The inclusion of the 
neglected terms lifts the degeneracy and gives a spectrum of two doublets u^2 and 
0^4. The corresponding amplitudes a a are free constants in a general solution which 
are determined by the initial conditions (only the ratios n a = b a /a a are fixed, as in 
(pD|)). It will turn out in the next section that a 3j4 are generally so small that co> 3j4 
cannot be observed. Only for very special conditions can W34 be seen. 

Finally we shortly discuss non-gyrotropic excitations which we define as having 
only even-order terms in the equations of motion. Examples are the kinks in the 
1-d nonlinear Klein-Gordon models for which G and A vanish. Here a 4 th -order 



equation was derived but not considered in detail [p]] . Another example are the 
planar vortices of the 2-d anisotropic Heisenberg model (stable for 5 > 5 C ). Here 
G = in Eq. (0) because of p = 0; all other odd-order terms also vanish because 
they are proportional to p, e.g. (|31~1). 



5 Comparison with simulations 

A single vortex with polar coordinates (R, 0) on a circle with radius L has an image 
vortex at Ri = L 2 /R and the same <fi (as in 2-d electrostatics). For free boundaries 
the image has opposite vorticity but the same polarization || and the force is F(R) = 
27r /(Ri — R). The equation of motion (|10|) has a steady-state solution: a constant 
rotation on a circle R = Rq, <p(t) = u>ot with the following relation between Rq and 

-Auj%Ro-Muj 2 Ro + Glu R = F(R ). (44) 

This circle is identified with the mean trajectory around which we observe oscillations 
in the simulations (details are given in the Appendix). In Fig. |]the radial coordinate 
R(t) is plotted vs. the coordinate Ro<f)(t) in azimuthal direction. This figure already 
shows qualitatively that two frequencies are involved. In fact, the Fourier spectra 
for r(t) = R(t) — Rq and (p(t) = <j)(t) — u> t in Fig. [5] clearly show two dominant 
frequencies 10^%. The phase shifts 8\^ are approximately tt/2 and —tt/2 (Table |l]). 
Thus the simulation data are described very accurately by 

r{t) = a x cos u>it + a 2 cos u; 2 t (45a) 

(45b) 



r{i) = ai cos u>ii ■+■ a 2 cos u; 2 t 
Ro<p(t) = bi sin uj\t + 6 2 sin co^ 



with 

signKo, = ±l , a = 1,2. (46) 
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As there are displacements not only in radial direction but also in azimuthal direction 
the trajectories are a superposition of two cycloids with frequencies uj\^ 2 . 

In order to compare with our theory we must solve the equation of motion (|T0D 
in polar coordinates. The result agrees completely with (|45), (|46|). The formula for 



uj a is much more complicated than the cartesian result fl3"8|) due to many additional, 
Wo-dependent terms. However, K a can be written in a rather simple form by using 



Kl = 1+ (M + 3 loK (47) 

f F> F ° - 2 * 2rf (AH) 

f ~ F °-R- o -V0^fy (48) 

with i] = Rq/L. For 77 «^C 1, which is the case in our simulations, a short calculation 
shows that n a = ±1, neglecting terms of order rf and higher. In this case it turns 
out that the substitution 

uj a ± u — > u a , a = 1, 2 (49) 

transforms the complicated eigenvalue equations for u a into simpler ones 

{G-Av 2 a )v a K a = Mvl + F^ (50a) 

(G-Au 2 a )u a = (Mv 2 a + F /R )K a . (50b) 

Here the difference between Fq and F /R must be neglected to be consistent with 



the above neglection of / in fl47l). We note in passing that this difference and the 
one between \ki\ and |«2| are responsible for the deviations of M and A from the 
symmetries in ( |32a| ) and (|32b| ), which are discussed below these equations. 



As we know only the size dependence of M and A : it does not make sense to solve 
)0a|) , (|50b| ) for v a . Instead, we calculate M and A as functions of z/ a , or u a and lj , 



which have been measured in simulations for different system sizes: 



with 



VI = 


(B 1 ul-B 2 u 3 1 )/D 


A = 


{B 2 k x v\ - B x k 2 vI)ID 


B a 


= Gv a — K a F Q / Rq 


D 


= v\v\{kxV 2 - K 2 Vi). 



(51a) 
(51b) 



(52a) 
(52b) 



The results in Table have to be compared with fl33|) , which was calculated for 
Ro = 0, however. Therefore we take only the data for small Ro/L = rj, e.g. those for 
1] w 0.17. The data for A are well represented by A = Aq + CL a with a = 2.002, 
C = 4.67 and A = 40. This agrees perfectly with the /^-dependence in (^). From 
C we obtain b = 2.97, which agrees rather well with the value 2.5 which was estimated 
at the end of section 3. However, the values for M are practically independent of 
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L, in contrast to the logarithmic dependence in (|33|). This can be explained by the 
following: In contrast to the l/r'-decay in (|26|), rri\ seems to go to an L-dependent 
constant at the boundary (however, this is difficult to see because only mi + m 2 
is observed, and m 2 is never exactly zero). At this point we must realize that Eq. 
(p6|) is only an approximation, because for O i n ( p5|) we have taken (0), i.e., the 



contribution of the image vortex has been neglected. This effect and others, such as 
the splitting of the Ma- components for Rq ^ and the influence of a second vortex 
in the system, will be considered in forthcoming work. 

The last point in this paper is the observation of additional frequency doublets 
which can be related to higher levels of the hierarchy (see the penultimate paragraph 
of section 4). We have performed simulations for a vortex in the center of a square 
system with antiperiodic boundary conditions. In this particular case there are many 
image vortices at positions which can be calculated as in electrostatics. The forces 
from all the images cancel exactly, therefore the vortex in the center is influenced only 
by the small pinning forces from the lattice. We observe that the vortex center moves 
on small-amplitude cycloidal trajectories, which are either completely inside a lattice 
cell or which go over severall cells. However, the spectrum is always practically 
the same which means that the pinning forces do not substantially influence the 
frequencies. In Fig. | we clearly see not only a second doublet but also a third one, 
which we want to relate to the 5 th and 7 th levels of the hierarchy, resp. In order 
to show this relation we must check that the parameters M and A do not change 
dramatically when we go from the 3 rd to the 5 th level, for instance. On this level we 
get a 4 th -order eigenvalue equation for ui a , after splitting off the trivial solution. We 
need not solve this equation because we want to calculate M, A, ... as functions of 
the u>aS. This is achieved by using Vieta's rules which express the coefficients of an 
n th -order polynominal equation by its roots. Table |3| contains the observed o; Q 's and 
the resulting values for M and A, which in fact do not differ much on the different 
levels of the hierarchy. We remark that this test is very sensitive: Small changes in 
u a result in large changes in the values of M and A. 

Finally we would like to stress that the additional uj^s can usually not be observed 
in the spectra because their amplitudes are too weak, thus our 3 rd -order equation of 
motion is normally sufficient to describe the simulations. 

6 Discussion 

Inserting the usual travelling-wave ansatz into the Hamilton equations shows that 

the shape of a nonlinear coherent excitation generally depends also on its velocity X. 

Iterating this process we have shown that there is an additional dependence on X, 

X, and so on. This yields a hierarchy of equations of motion for X(t). For the case of 
the non-planar vortices of the 2-d anisotropic Heisenberg model the odd levels of this 
hierarchy represent an increasingly better description of the dynamics as observed in 
computer simulations of the full spin system. 

13 



Naturally the question arises why spin waves do not appear here, neither in theory 
nor in simulations. As to the latter point, we make two remarks: 1.) For a short time 
after the start of a simulation spin waves are radiated because the initial condition 
usually does not perfectly represent a moving vortex. In fact, if we use approximate 
formulas for the vortex structure, the generated spin waves are clearly observed. In 
our older simulations 0, || we eliminated them by adding a Gilbert damping term 
to the Landau-Lifshitz Eq. for the above short time period. However, in the present 
paper we use an iterative method (see Appendix) which produces a stationary vortex 
solution which is so good that practically no spin waves are generated. 2.) One might 
think that a vortex would continously radiate spin waves because it is subject to 
accelerations on its complicated cycloidal trajectory. Even if the amplitudes of these 
spin waves were too small to observe them directly, the effect should be detected 
indirectly by an energy loss of the vortex. Interestingly, in our simulations, even 
after a long time (several periods T = 27i/u) , see (ffl|)) there was no detectable 
energy loss; this suggests the possibility of additional conservation laws. 

As to the theory, there is an alternative formalism |T2[ which includes spin waves. 
Instead of eqs. (|8a|), (|8b|) the following ansatz is used 

0(f,t) = <j>(f-Y(t))+£(f,t) (53a) 

m(f,t) = m(f—Y(t))+x(r,t) (53b) 

where the functions <fi and m on the r.h.s. represent the static structure of a vortex, 
while £ and x represent spin waves. As the number of degrees of freedom on the 
r.h.s. of eqs. (^) is larger than on the l.h.s., constraints must be introduced which 



have the form of orthogonality relations [13, 11 



One of the constraints leads to an implicite definition of the vortex position Y. 
This definition is global in the sense that all spins in the system are involved, in 
contrast to the local definition of the vortex center X that has been used in the 
present paper (see Appendix). Consequently different trajectories Y(t) and X(t) are 
obtained by analyzing the same simulation data. Y(t) turns out to be equivalent to 
the mean trajectory Ro(t) belonging to the cycloidal trajectories X(t) obtained in 
the present paper. 

In the alternative formalism the cycloidal trajectories are obtained by a coupling 
of Y(t) to certain quasi-local spin wave modes. Far away from the vortex center 
these modes are extended while they are localized in the core region 14]. Quasi- 



local eigenmodes are in contrast to the truly localized, intrinsic modes of e.g. the 
kinks in the 4 -model [|15"1| . An exact numerical diagonalization for a small system 



(L = 20) with one vortex shows that essentially only two modes are involved. Their 
frequencies are identical to those of the doublet which we observe in the spectrum 
of the cycloidal trajectories (e.g., Table [I] and Fig. ||). If higher-order doublets can 
be seen in the spectrum (like in Fig. |]), they can also be identified with spin waves. 
Hence, going e.g. from the 3 rd to the 5 th level of the hierarchy in the present paper 
corresponds to taking into account two additional spin wave modes in the alternative 

14 



formalism. This strong model selection may also be responsible for inhibiting spin 
wave emission by the vortex motion as noted above. 

Acknowledgements 

Work at Los Alamos National Laboratory was sponsored by the United States De- 
partment of Energy. Travel between Los Alamos and Bayreuth was partially sup- 
ported by NATO Collaborative Research Grant No. 0013/89. 

Appendix 

In this appendix we describe the numerical details of our work. 

The time evolution of classical spin systems is governed by the Landau-Lifshitz 

equation 

a tt 

S k = S k xH k , H k = - (54) 

dS k 

The index k denotes the site-coordinates of a 2-d square lattice. The local field H k 
is for our model (§) given by 

Hk = JYXSfZ* + S ?ey + C 1 - OT^] ( 55 ) 

(ki) 

where here the sum runs over all the nearest neighbors of k not including k itself. 

Equation fl54|) with (|55|) was numerically integrated using a forth-order Runge- 
Kutta scheme with a time step of 0.04 (in time units (JS 1 ) -1 ). Either free bound- 
ary conditions on a circular system, i.e. Sij = for all (i,j) with i 2 + j 2 > L 2 , 
or anti-periodic boundary conditions on a square system were used. With 
tan -1 Sfj/Sfj the latter are defined in the following way: 



'i,i 



0i,O + <f>i,L = V0 , <Pi,l + <Pi,L+l = Vo 

4>0,j + 4>L,j = ¥0 + QinTT , 4>l,j + 4>L+l,j = V?0 + qin^ 

$1,0 + $i,L = ) <Jj,l + Si,L+l = 

0,j + ^L,j = > "l,j + ^L+l,j = 0- 

L denotes the number of lattice sites in x- and ^/-direction, respectively (i.e. the linear 
dimension of the system). q in is the total amount of vorticity in the system (1 or —1 
for one vortex in the system) and (fo is an arbitrary constant. 

The exact structure of an one-vortex solution is, even in the static case, analyt- 
ically not known. Therefore we calculated the initial conditions for our simulations 
numerically using the following iteration scheme: 

I c?( n ) I 
Sf +1 > = ^-H^. (56) 

\ n k 
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Normally this method is used to calculate static solutions of the Landau-Lifshitz 
equation. However, iterating only a few times (a typical value is 30) a pronounced 
vortex structure develops which is suitable to serve as the initial spin field for the 
simulations. 

Next we explain how we have determined the vortex position(s) during a simula- 
tion. Using the discrete variant of the contour integral §drV(f)(r) (which yields 27rg if 
the contour surrounds the vortex and zero otherwise), only the lattice plaquette can 
be identified in which the vortex necessarily has to reside. We estimated the precise 
position according to the formulas 

v — I - sinfc-^Qi+^-tfe-tM)] (57a") 



2 cos[g- 1 (</)i+(/)2-(/<3-'/>4)]-sin[g- 1 ((?ii-fli2+'?i3-</ , 4)]-cos[g- 1 ((/)i-(/)2+</i3-04)] 

X _ 1 sin[iji' :L (0i-02-03+04)] f57b) 

2 2 COs[(J- 1 ((/)i-02-</'3+</'4)]+sin[(?~ 1 (?il-?i2+</'3-04)]- cos [9" 1 ('/ , l-</ , 2+03-</'4)] ^ ' 

where 1; . . . , 4 label the angles (cf. (|T5|a)) of the four innermost spins at the vortex 
core, beginning with the spin downleft and surrounding the vortex counter-clockwise. 
X\ and X 2 are measured in units of the lattice constant. The underlying coordinate 
system has its origin at the center of the lattice plaquette under consideration, i.e. 
X\ and X 2 range between [—0.5, 0.5] for 'reasonable' values of (f>i, . . . , 04. 

Our basic assumption in deriving fl57|) was that 0i, . . . , 04 are distributed accord- 
ing to the static solution ([18]). Considering only differences, this eliminates the rather 
annoying constant <y9 in (|18"D, we can derive the following four equations: 



i i _i_ oy 
tan[g- 1 (0 1 -0 2 )] = -- 2 (58a) 



2 xi + Xi + X 



2 



1 1 - 2Xt 

2 J\^ -\- X2 — X\ 



tan[g- 1 (0 2 -0 3 )] = ~ v , v2 X — (58b) 



-\(± m l 1 ~ 2X 2 



tan[g (0 3 - 4 )] = ~2 X 2 + X 2 -X ^ 8 ^ 



tan[g- 1 (0 4 -0 1 )] = -- v2 - y2 I v . (58d) 



1 1 + 2Xi 

2xJTxJTx 1 

We formally assume that the quantities X\, X 2 and Xf + X\ are independent of 
each other. Eqs. (^b,d) then constitute a linear system in X\ and X\ + Xf and 
analogeously (|58|a,c) in X 2 and X\ + X|. Both set of equations can easily be solved 
and yield after some further trigonometric steps the above formulas for X\ and X 2 . 
Of course in general, 0i, . . . , 04 are not distributed according to the tan - ^function. 
Although then (|58|a-d) do (most probably) not have a solution at all, the Eqs. (|57|) 
can still be justified in the sense that they provide an approximation of the minimal 
residual of (|58"|a-d) . Anyway, the reasonability of (^7|) is also acknowledged by the 
numerical simulations which show that trajectories are smoothly connected when the 
vortex leaves one lattice plaquette and enters another one. 
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1.04 - 




Figure 1: Out-of-plane structure of the vortex at the 7 th turning point of the trajec- 
tory in Fig. [|. Here the acceleration has a maximum and points in radial direction, 
while the velocity is small and points in azimuthal direction. 




Figure 2: Out-of-plane structure at the middle between the 7 th and the 8 th turning 
point of Fig. [|. Here the velocity has a maximun and points in radial direction, while 
the acceleration has a minimum. 
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Figure 3: Same as in Fig. p], but at the 8 th turning point, where the accelaration 
points in negativ radial direction. 




Figure 4: First part of the trajectory of a vortex with q = p = 1 on a circular system 
with radius L = 36 and free boundary conditions. The small diamonds (o) mark the 
position of the vortex in time intervals of 10( JS)^ 1 . 
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Figure 5: (a) Fourier spectrum of the radial displacements r(t) = R(t) — Rq, from 
simulation data for < t < 20000(JS')~ 1 . The first part of the trajectory is plottet 
in Fig. 01. (b) Spectrum of the aximuthal displacements tp(t) = <p(t) — uJot. 
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Figure 6: Spectrum of one of the cartesian coordinates of a vortex on a 50 x 50 
square lattice with antiperiodic boundary conditions, from simulation data for < 
t < 20000(JS)- 1 . 
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Table 1: Observed data from the trajectory of a vortex with q = p = 1 on a circle of 

L R cj /10- 3 cji/10- 2 w 2 /10- 2 <*i <5 2 «i k 2 



24 


4.03 


1.84 


4.257 


5.302 


1.57 


-1.57 


1.01 


-1.06 


24 


6.02 


1.90 


4.257 


5.276 


1.56 


-1.56 


1.03 


-1.11 


24 


8.02 


2.01 


4.261 


5.249 


1.56 


-1.56 


1.06 


-1.18 



36 


5.80 


0.814 


2.958 


3.426 


1.57 


-1.57 


1.00 


-1.02 


36 


8.28 


0.823 


2.962 


3.417 


1.52 


-1.50 


1.02 


-1.09 


36 


9.87 


0.852 


2.964 


3.412 


1.56 


-1.56 


1.06 


-1.12 


36 


12.25 


0.889 


2.966 


3.404 


1.56 


-1.57 


1.05 


-1.14 


72 


16.11 


0.205 


1.546 


1.661 


1.43 


-1.42 


- 


- 


72 


24.20 


0.218 


1.548 


1.657 


1.39 


-1.40 


1.12 


-1.22 


72 


31.95 


0.243 


1.550 


1.653 


1.45 


-1.46 


1.24 


-1.37 



Table 2: Parameters M and A of the equation of motion, calculated from ( |51a| ), 
( |51bD . Values in parentheses are extrapolated. 

L R r] = R /L M A 



24 


4.03 


0.168 


13.7 


2750 


24 


6.02 


0.251 


12.5 


2764 


24 


8.02 


0.334 


10.8 


2776 



36 


5.80 


0.161 


13.9 


6167 


36 


8.28 


0.230 


12.9 


6174 


36 


9.87 


0.274 


12.0 


6180 


36 


12.25 


0.340 


10.7 


6190 


72 


12.0 


0.167 


(13.8) 


(24416) 


72 


16.1 


0.224 


13.0 


24436 


72 


24.2 


0.336 


10.7 


24477 


72 


31.9 


0.443 


7.4 


24504 
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Table 3: Parameters M and A of the third-order equation of motion, compared to 
the parameters M, A, B, C of the fifth-order equation, where B and C belong to 
X^> and X*- 5 -*, respectively. The cj q 's are observed in vortex motion on 2L x 2L 
square lattices with antiperiodic boundary conditions. 

L 10 25 



cji/10" 1 


1.376 


0.5435 


c^/lCT 1 


1.696 


0.5938 


CO3/IO- 1 


3.749 


1.2692 


W4/IO- 1 


4.362 


1.3670 


M (3 rd ) 


8.62 


9.79 


A (3 rd ) 


269 


1947 


M (5 th ) 


10.98 


13.42 


A (5 th ) 


304 


2300 


B (5 th ) 


153 


1689 


C (5 th ) 


1646 


111991 
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